function [y, bb] = aft_SampleBox(hfile, bb, vol, imeth, trans)
% AFT::SampleBox  - sample VoxelData within a BV based box
%
% FORMAT:       y = obj.SampleBox([bb [, vol [, imeth [, trans]]]])
%
% Input fields:
%
%       bb          BoundingBox struct with fields
%        .BBox      XYZ start/end, default is a MNI "default" box
%        .XYZRes    resolution (either 1x1 or 1x3), default: [3, 3, 3]
%       vol         volume number (default 1)
%       imeth       method {'linear'}, 'cubic', 'nearest'
%       trans       optional 4x4 transformation matrix applied to the
%                   coordinates
%
% Output fields:
%
%       y           double precision data sampled at given coordinates
%
% Note: This method requires the MEX file flexinterpn. Also, for use
%       in BVQX (i.e. as VTC data), the resolution must be uniform!

% Version:  v0.8a
% Build:    9102722
% Date:     Oct-27 2009, 10:15 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
sc = bvqxfile_getscont(hfile.L);
st = sc.S.Extensions{1};
if ~any(strcmpi(st, {'ava', 'glm', 'hdr', 'head', 'vdw', 'vmp', 'vmr', 'vtc'}))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid file type for ::SampleBox.' ...
    );
end
if nargin < 2 || ...
   ~isstruct(bb) || ...
    numel(bb) ~= 1
    bb = struct;
end
if ~isfield(bb, 'BBox') || ...
   ~isequal(size(bb.BBox), [2, 3]) || ...
   ~isa(bb.BBox, 'double') || ...
    any(isinf(bb.BBox(:)) | isnan(bb.BBox(:)) | bb.BBox(:) < 0 | bb.BBox(:) > 256) || ...
    any(diff(bb.BBox) < 0)
    bb.BBox = [44, 38, 44; 241, 193, 211];
end
if ~isfield(bb, 'XYZRes') || ...
   ~isa(bb.XYZRes, 'double') || ...
   ~any([1, 3] == numel(bb.XYZRes)) || ...
    any(isinf(bb.XYZRes) | isnan(bb.XYZRes) | bb.XYZRes < 1 | bb.XYZRes > 15 | bb.XYZRes ~= fix(bb.XYZRes))
    bb.XYZRes = [3, 3, 3];
else
    if numel(bb.XYZRes) ~= 3
        bb.XYZRes = bb.XYZRes([1, 1, 1]);
    end
end
if nargin < 3 || ...
    isempty(vol)
    vol = 1;
elseif numel(vol) ~= 1 || ...
   ~isa(vol, 'double') || ...
    isinf(vol) || ...
    isnan(vol) || ...
    vol < 1 || ...
    vol ~= fix(vol)
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid vol argument.' ...
    );
end
if nargin < 5 || ...
    isempty(trans)
    trans = eye(4);
elseif ~isa(trans, 'double') || ...
   ~isequal(size(trans), [4, 4]) || ...
    any(isinf(trans(:)) | isnan(trans(:))) || ...
    any(trans(4, :) ~= [0, 0, 0, 1])
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid trans argument.' ...
    );
else
    trans = [0, 0, 1, 0; 1, 0, 0, 0; 0, 1, 0, 0; 0, 0, 0, 1] * trans * ...
        [0, 1, 0, 0; 0, 0, 1, 0; 1, 0, 0, 0; 0, 0, 0, 1];
end

% for non-BV types, get coordinate frame, otherwise the bounding box
if strcmpi(st, 'head')
    cfr = head_CoordinateFrame(hfile, vol);
elseif strcmpi(st, 'hdr')
    cfr = hdr_CoordinateFrame(hfile, vol);
else
    mbx = aft_BoundingBox(hfile);
    cfr = struct('Trf', mbx.QuatB2T);
end
try
    tmat = inv(double(cfr.Trf)) * trans;
catch
    error( ...
        'BVQXfile:InvalidArgument', ...
        'Coordinate transformation matrix not invertible.' ...
    );
end
if nargin < 4 || ...
   ~ischar(imeth) || ...
   ~any(strcmpi(imeth(:)', {'cubic', 'lanczos3', 'lanczos8', 'linear', 'nearest'}))
    imeth = 'linear';
else
    imeth = lower(imeth(:)');
end
bb.BBox(2, :) = bb.BBox(1, :) + bb.XYZRes .* ceil(diff(bb.BBox) ./ bb.XYZRes) - bb.XYZRes;
bb.SBox = [ ...
    Inf, Inf, Inf; ...
    128 - bb.BBox(1, [3, 1, 2]); ...
    -bb.XYZRes; ...
    128 - bb.BBox(2, [3, 1, 2])];

% get data in single precision
try
    vdt = single(aft_GetVolume(hfile, vol));
catch
    rethrow(lasterror);
end

% set NaNs to 0 first
vdt(isnan(vdt)) = 0;

% sampling
tmatb = tmat([2, 3, 4, 5, 7, 8, 9, 10, 12]);
if all(abs(tmatb) <= eps)
    bb.SBox(2:4, :) = bb.SBox(2:4, :) .* tmat([1, 6, 11; 1, 6, 11; 1, 6, 11]);
    bb.SBox([2, 4], :) = bb.SBox([2, 4], :) + tmat([13, 14, 15; 13, 14, 15]);
    tmat = {};
    bb.Trf = eye(4);
else
    bb.Trf = tmat;
    tmat = {tmat};
end
y = permute(flexinterpn_method(vdt, bb.SBox, 0, tmat{:}, imeth), [2, 3, 1]);
